
function paiOpt = GetTimingPrem(outEstim,pai0,numSim,pathLength)

model  = outEstim.model;
params = model.params;

% For anti-thetic sampling
rng(1,'twister');
shocksMat = zeros(size(model.eta,2),pathLength,numSim);
shocksMat(:,:,1:numSim/2) = randn(size(model.eta,2),pathLength,numSim/2);
shocksMat(:,:,numSim/2+1:end) = -shocksMat(:,:,1:numSim/2);

% settings
setupTP.model  = model;
setupTP.params = model.params;
posVf          = find(strcmp(model.labely,'$mvf_t$'));
setupTP.vf     = exp(model.g0(posVf) + 1/2*model.gss(posVf) + 1/6*model.gsss(posVf));


% simulate
posC      = find(strcmp(model.labely,'$c_t$'));
posW      = find(strcmp(model.labely,'$w_t$'));
pos_muz   = find(strcmp(model.labelx,'$\mu _{z,t}$'));
pos_n     = find(strcmp(model.labelx,'$n_{t}$'));
pos_d     = find(strcmp(model.labelx,'$d_{t}$'));

simC = NaN(numSim,pathLength);
simZ = NaN(numSim,pathLength);
simD = NaN(numSim,pathLength);
simN = NaN(numSim,pathLength);
simW = NaN(numSim,pathLength);

for i = 1:numSim
    %shocks = randn(size(model.eta,2),pathLength);
    shocks = shocksMat(:,:,i);
    [ySim,xfSim,xsSim,xrdSim] = simulate_3rd_kron_v2(model.gx,model.gxx,model.gxxx,...
                                                         model.hx,model.hxx,model.hxxx,...
                                                         model.gss,model.gssx,model.gsss,...
                                                         model.hss,model.hssx,model.hsss,...
                                                         model.eta,shocks,outEstim.setupFilter.appOrder);    
    xSim = xfSim + xsSim + xrdSim;
    
    muz       = model.h0(pos_muz)   + xSim(pos_muz,:);
    n         = model.h0(pos_n)     + xSim(pos_n,:);
    d         = model.h0(pos_d)     + xSim(pos_d,:);
    simZ(i,:) = cumprod(exp(muz),2); 
    simN(i,:) = exp(n);
    simD(i,:) = exp(d);
    simC(i,:) = exp(model.g0(posC) + ySim(posC,:)).*simZ(i,:);
    simW(i,:) = exp(model.g0(posW) + ySim(posW,:)).*simZ(i,:);
end

bettap = params.BETTA.^(0:1:pathLength-1);

setupTP.N0  = params.U0*sum(repmat(bettap(1,1:end-1),numSim,1)...
              .*simZ(:,2:end).^((1-params.CHI)*(1-params.CHI0)),2);
Css         = exp(model.g0(posC));
setupTP.Nc  = 1/(1-params.CHI)*sum(repmat(bettap(1,1:end-1),numSim,1).*simD(:,2:end)...
              .*((simC(:,2:end)-params.B*simC(:,1:end-1))./((Css*simZ(:,2:end)).^params.CHI0)).^(1-params.CHI),2);

setupTP.Nd  = params.U0d*sum(repmat(bettap(1,1:end-1),numSim,1).*simD(:,2:end)...
              .*simZ(:,2:end).^((1-params.CHI)*(1-params.CHI0)),2);

PHIzero     = model.params.PHIzero;
setupTP.Nl  = PHIzero^(params.PHI)/(1-1/params.PHI)*...
              sum(repmat(bettap(1,1:end-1),numSim,1).*simD(:,2:end).*simN(:,2:end).^params.PHI...
             .*simZ(:,2:end).^(params.PHI*(1-params.CHI)*(1-params.CHI0))...
             .*(simC(:,2:end)-params.B*simC(:,1:end-1)).^(params.CHI*(params.PHI-1))...
             .*(Css.*simZ(:,2:end)).^(params.CHI0*(params.PHI-1)*(1-params.CHI))...
             ./simW(:,2:end).^(params.PHI-1),2);
setupTP.c_cu           = Css;
setupTP.w_cu           = exp(model.g0(posW));
setupTP.z_cu           = 1;
setupTP.n_cu           = 1;
setupTP.d_cu           = 1;
setupTP.params         = params;
setupTP.params.PHIzero = PHIzero;
setupTP.Css  = Css;
% Q = timingPremiumObjectFunc(0.5,setupTP);

options = optimset('Display','iter','MaxFunEvals',2000,'MaxIter',2000,'TolFun',1D-4,'TolX',1D-3);
paiOpt = fminsearch(@timingPremiumObjectFunc,pai0,options,setupTP);
%paiOpt    = fminunc(@timingPremiumObjectFunc,pai0,options,setupTP);

end